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Abstract 

The (111), (HO), and (001) surfaces properties of PuC>2 are studied by using density-functional 
theory+C/ method. The total-energy static calculations determine the relative order of stability 
for low-index PuC>2 surfaces, namely, O-terminated (111) > (110) > defective (001) > polar (001). 
The effect of thickness is shown to modestly modulate the surface stability and chemical activity of 
the (110) surface. The high work function <£> of 6.19 eV indicates the chemical inertia of the most 
stable (111) surface, and the surface O-vacancy with concentration Cy = 25% can efficiently lower 
<£> to 4.35 eV, which is a crucial indicator of the difference in the surface chemical activities between 
PuC>2 and a-Pu203. For the polar (001) surface, 50% on-surface O-vacancy can effectively quench 
the dipole moment and stabilize the surface structure, where the residual surface oxygen atoms are 
arranged in a zigzag manner along the <100> direction. We also investigate the relative stability 
of Pu02 surfaces in an oxygen environment. Under oxygen-rich conditions, the stoichiometric O- 
terminated (111) is found to be the most stable surface. Whereas under O-reducing conditions, 
the on-surface O-vacancy of Cy = 1/9 is stable, and for high reducing conditions, the (111) surface 
with nearly one monolayer subsurface oxygen removed (Cy = 8/9) becomes most stable. 

PACS numbers: 71.15.Mb, 71.30. +h, 71.28.+d, 71.27.+a 
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I. INTRODUCTION 



Plutonium dioxide (PUO2) is of high importance in the nuclear fuel cycle and is par- 
ticularly crucial in long-term storage of Pu-based radioactive waste. Besides playing an 
important role in both technological applications and environmental issues, Pu metal and 
its oxides show many intricate physical behaviors due to the complex electronic structure 
properties of 5f states jl-4|. Therefore, a thorough understanding of the physical and chem- 
ical properties of PuC>2 is of great significance, full of challenges and has always attracted 
attention. Recently, there have occurred in the literature a series of experimental reports 
5KZJ on the strategies of storage of Pu-based waste. When exposed to air and moisture, 
metallic plutonium surface rapidly oxidizes to PuC>2 Under special aqueous condi- 

tion, the interaction of PuC>2 surface with adsorbed water can generate non-stoichiometric 
Pu0 2 + x (x<0.27) {5] via an overall reaction, namely, PuC^+xH^O— >Pu02+ x +H 2 . However, 
the oxidation of PuC>2 has been proved to be strongly endothermic by subsequent first- 
principles theoretical calculation |10l ? ]. To test the possible existence of surface PuC^+x, 
recent photoemission study [7] has been carried out and found that Pu0 2 is only covered by 
a chemisorbed layer of oxygen and can be easily desorbed at elevated temperature. Thus, 
PuC>2 is generally acknowledged as the highest stable Pu-oxide under ambient conditions. 
However, under oxygen- lean conditions (in the vacuum or inert gas), the PuCVlayer can 
be reduced to sesquioxides (Pu 2 3 ), which can promote the corrosion of the Pu-metal by 
hydrogen jsj]. As we know, the low-temperature phase of Pu-sesequioxide is a phase with 
space group Ia3 (No. 206), which is similar in the crystal structure to the cubic PuC>2 
(Fm3m, No. 225) with the 25% O vacancy located in the 16c (0.25, 0.25, 0.25) sites. The 
above mentioned experimental observations indicate that the surface of PuC>2 is to some 
extent chemically inactive, however, the formation of O vacancies can prominently modify 
the electronic structure properties of both the bulk and the surface of PuCV As a matter 
of fact, the surface layers are directly involved in the significant corrosion processes and 
many technological applications of the actinide oxides, thus a deep understanding of the 
physical and chemical properties of PuC>2 surfaces is always desirable. However, due to the 
radioactivity and toxicity of Pu and the complexity of the Pu element and Pu-0 system, 
it is extraordinarily difficult to experimentally explore the surface atomic and electronic 
structure properties of the Pu-oxides, and particularly so for a single phase compound. 
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From the theoretical point of view, conventional density-functional theory (DFT) that 
applies the local density approximation (LDA) or generalized gradient approximation (GGA) 
underestimates the strong on-site Coulomb repulsion of the 5f electrons and, consequently, 



describes PuO"2 as incorrect ferromagnetic FM conductor instead of antiferromagnetic 
AFM Mott insulator reported by experiment [121 ] . Similar problems have been confirmed 
in studying other correlated materials within the pure LDA/GGA schemes. Fortunately, 



several approaches, including the LDA/GGA+c/ the hybrid density functional of 

(Heyd, Scuseria, and Enzerhof) HSE [161 ]. the self-interaction corrected local spin density 



SIC-LSD |17|, and LDA combined dynamical mean- field theory DMFT 18], have been 
developed to correct these failures in calculations of actinide compounds. The effective 
modification of pure DFT by LDA/GGA+f/ formalisms has been confirmed widely in study 



of Pu0 2 



19 



31-42-41] . By tuning the effective Hubbard parameter in a reasonable range, the 
AFM Mott insulator feature was correctly calculated and the structural parameters as well 
as the electronic structure are well in accord with experiments. However, those increasing 
theoretical researches have been focusing on the bulk properties of Pu02, and very little 
is known regarding its surface physical and chemical properties, which is in sharp contrast 
to the depth and comprehensiveness of researches conducted upon the transition metal and 

^T studies of the Pu0 2 (100) and 



In addition, one of the most 



rare earth oxides [25(. As far as we are aware, few D 
(110) surfaces have been presented in the literature {26 . 
important issues, i.e., the possible formation of O vacancies and their effect on the atomic 
and electronic structures of Pu0 2 surfaces, remains completely unexplored. 

Motivated by the above mentioned facts, in this paper, we systematically study the surface 
properties of Pu0 2 . Specifically, we have addressed (i) the structural stabilities of the low- 
index Pu0 2 (lll), (HO), and (001) surfaces based on the calculations of surface energies and 
surface relaxations, (ii) the surface electronic properties such as the surface work function 
and layer-projected density of the state (PDOS), and (iii) the effects of Pu-oxide thickness 
and O-vacancy on the surface stabilities and electronic structure properties. 

This paper is organized as follows. The details of our calculations are described in Sec. 
II. In Sec. Ill we present and discuss the results. In Sec. IV, we summarize our main 
conclusions. 
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II. METHODOLOGY OF THE CALCULATION 



DFT calculations 



The DFT calculations are carried out using the Vienna ab initio simulation 



28[ with the frozen-core projector-augmented wave (PAW) pseudopotentials |29|, |30| and 



package 



plane- wave set. For the plane- wave set, a cut-off energy of 400 eV is used. The plutonium 
6s 2 7s 2 6p 6 6d 2 5f 4 , and the oxygen 2s 2 2p A electrons are treated as valence electrons. The ex- 
change and correlation effects are treated in both the LDA and the GGA 31], based on 
which the strong on-site Coulomb repulsion among the localized Pu 5f electrons is described 
by using the LDA/GGA+[7 formalisms formulated by Dudarev et a/.[13Nl5|. As concluded 



in some previous DFT studies [19|, 



20 



22 



24] . although the pure LDA and GGA fail to 



depict the electronic structure, especially the insulating nature and the occupied-state char- 
acter of bulk AFM Pu02, the LDA/GGA+t/ approaches can capture the Mott insulating 
properties of the strongly correlated Pu 5/ electrons adequately and well reproduce experi- 
mental ground-state parameters by tuning the effective Hubbard parameter U e g at ~4 eV. 
In this paper, the spherically averaged screened Coulomb energy U and the exchange energy 



J for the Pu 5/ orbitals are set to be 4.75 and 0.75 eV respective 



y, which have been tested 



and applied in our previous studies of plutonium oxides [19|, |20|, |24| . For fluorite structure 
Pu02 of AFM, our calculated equilibrium lattice parameter ao=5.466 A within GGA+£7 or 
an=5.362 A within LDA+f/ is in good agreement with the experimental value of 5.398 A 
[5] . Our extensive test calculations in this work indicate that the choice of U e g can alter the 
electronic-structure properties of Pu0 2 surfaces, as well as those of bulk Pu0 2 . Specifically, 
when U e ff is less than 2.0 eV, the results of the electronic density of state (DOS) indicate 
that PUO2 surfaces are metallic FM-conductor instead of the AFM Mott-insulators. The 
combination of 11=4.75 and J=0.75 eV is also the optimum to well describe the electronic 
structure properties of the surfaces of PUO2, although the atomic structural optimizations 
seem to be insensitive to the choice of U e g. 

The low-index surfaces of PUO2 are modeled by finite-sized periodic supercells, consisting 
of a number of oxide layers infinite in x and y directions and separated in the z direction by 
a vacuum of 30 A. The Brillouin zone (BZ) integration is performed using the Monkhorst- 
Pack (MP) Appoint mesh 32J. Specifically, the two-dimensional (2D) unit (1 x 1) cell with 
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(Ill) (110) (001) defect-(OOl) 




FIG. 1: (Color online) Low-index PuG"2 surface models: (a) ideal oxygen-terminated (111) unit 
cell; (b) (110) unit cell; (c) (y/2x ^)RA5° super-cell of ideal oxygen-terminated (001); and (d) the 
defective (001) surface with four different distributions of 50% O-vacancy on the bottom and top 
faces, namely, surf-1, surf-2, surf-3 and surf-4. The blue, red and white spheres denote Pu-atom, 
O-atom and O-vacancy, respectively. In the side-views of three ideal surface, the labels of different 
Pu-cation/O-anion [+ + /— ] stacking sequences are also indicated at the lower part of the graphs. 

a 11 x 11 x 1 Appoint MP grid is employed for the defect-free (ideal) PuC>2(lll) and (110) 
surfaces, and a (a/2 x y / 2)i?45° supercell with a 7 x 7 x 1 fc-point grid is used for the (001) 
surface. Examples of the initial slab configurations used in the calculation are shown in 
Figs. l(a)-(c), which are obtained by truncating bulk PuC>2 along the [111], [110] and [001] 
orientations, respectively. In the spin-polarized calculations with the AFM order set to be 
in a simple "t i t V alternative manner along the z direction, all atoms are fully relaxed 
until the Hellmann-Feynman forces are less than 0.01 eV/A. Various convergence tests have 
been performed to ensure the above mentioned input parameters and models feasible and 
reasonable in our current calculations. The result thereof shows that the surface energy of a 
slab with certain thickness is converged within 3 me V/A 2 . In this paper, for the convenience 
of depiction and plotting, we are using the number of Pu-cation layer N in a PuC>2 slab to 
represent its thickness. Thus, one can see in Figs. l(a)-(c) that the ideal (111), (HO), 
and (001) slabs have the same thickness tag, namely, N= 6, and it is worth noting that in 
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practice these slabs consist of 18, 6 and 13 atomic monolayers (MLs), respectively. 

According to Tasker's conclusion j^] on the surface stabilities of ionic crystals, the side- 
view of three low-index PuC>2 surfaces in Fig. 1 can reveal that the oxygen-terminated (111) 
surface, consisting of successive and electrically neutral "O-Pu-O" blocks, is the "Tasker- 
type-II " surface, the (110) is the typical "type- 1 " surface stacked with identical neutral 
planes fulfilling the PuC>2 stoichiometry, and the (001) is the typical "type-Ill " surface, 
consisting of oppositely charged planes. Generally, type-I and type-II are stable nonpolar 
surfaces, while type-Ill is unstable polar surface due to the dipole moment of the repeated 
unit in z direction. The dipole moment may be quenched and the polar surface stabilized 
by a variety of methods 34 3^] , including surface reconstruction, the presence of adsorbates, 
and changes in the surface electronic structure. In our present work, the polar (001) surface 
is modeled as two oxygen-terminated surfaces with 50% oxygen vacancies to fulfill the stoi- 
chiometric formula, and the different distributions of the 50% O-vacancy considered here are 
showed in Fig. 1(d). Besides above-mentioned three low- index surfaces, the Pu-terminated 
(111) surface generated by removing the outmost O layers of the O-terminated (111) sur- 
face has been found to be quite unstable and eventually become O-terminated through a 
significant reconstruction in the surface region due to the intense dipole-dipole interaction. 
In such stable structure, the upper layers resemble the /3-Pu2O3(0001) surface. Thus, we 
exclude the Pu-terminated (111) slab model and briefly name the O-terminated Pu02(lll) 
surface as (111) surface if not mentioned differently. 

The surface energy E s is the energy needed to produce a unit surface from a 3D infinite 
crystal and is one central quantity in the studies of the relative stability of different surfaces. 
In the DFT total-energy calculations of repeated slab-supercell geometries, E s can be written 
as 

77>relax/unrelax ^ ( /-.relax/unrelax „ \ /-. n 

= ^slab -^bulkj, (1) 

where ^^ ax / unrelax j s ^he total energy of the supercell with relaxed/unrelaxed slab, -Ebuik is 
the energy of the reference bulk with the same number of atoms, and the denominator 2A 
is the total area of both surfaces of the slab with a finite thickness. Here, the convergence 
of the E s (i.e., E™ in Eq.(l)) is mainly determined by two correlated factors, namely, 
the atomic structural relaxations in several outmost layers and the thickness of the slab 
model. For the surface structural relaxations, one can simply evaluate its contribution to 
E s by calculating the surface relaxation energy AE S = — (i^ elax — ^ nrclax ). The effect of 



the slab thickness should be highlighted especially when the slab consists of the complicated 
compounds such as the actinide oxides, and in addition, the work functions calculated with 
slab approximations are known to be depending on the slab thickness. Thus, in this work, 
both factors will be considered and discussed in detail. 



B. Thermodynamic considerations 

The DFT total-energy calculation gives E s only at zero temperature T=0, zero pressure 
P—0, and for the surface in contact with vacuum, which cannot be used to study the 
influence of the realistic environmental conditions at a specific T and P. To further study 
the relative stability of the PuC>2 surfaces with various concentrations of surface vacancy 
(Cy) at finite T and gas partial P of the surrounding environment, we take the approach 



of u ab initio atomic thermodynamics" 38|, |39| to get the surface Gibbs free energy 7(T,P), 
with the general expression given by 

1 



^ 2A 



G(T,p,{m})-J]T^iCz\pO 



(2) 



where G is the Gibbs free energy of the solid with the surface of interest, 2 A is the total 
surface and pi are the particle number, the chemical potentials and the partial 

pressures of the various species. Here, the focus of our work is the relative stability of O- 
terminated Pu02 surfaces with different O-vacancy concentrations, thus only two chemical 
species need to be considered, namely i=Pu and O. In practice, the surface Gibbs energy 
difference Aj(T, P) between a defective Pu02 surface and the corresponding defect-free 
(ideal) surface is the important quantity required, which can be written as 

A 7 (T, P) = i [G dcfcct (T, p, Ny) - G idcal (T, P) + N^ (T } p 02 )} , (3) 

where G defect and G ldeal are the Gibbs free energies of the supercells with the defective and 
ideal Pu02 surfaces, respectively, and Ny is the total number of O vacancies on the Pu02 
surface. In the present situation, the entropy and volume effects are small compared to the 
band energy in the Gibbs free energy and thus are neglected in our calculations. Ho{T,Po 2 ) 
in Eq. (3) is the oxygen chemical potential under partial pressure po 2 and for ideal oxygen- 
gase we can use the well-known thermodynamic expression 



/io(T,p 02 ) = \ (E O2 +fio 2 (T,p ) + k B T\n(p Jp )) , (4) 



where Eq 2 is the total energy of the oxygen molecule. For the standard pressure p° = latm, 
the values of flo 2 (T,p ) have been tabulated in Ref. 40. If we refer the fj,o to \Eq 2) then 
the allowed range for the Afx = /i — \E , 2 is given by 

-\e,^ A/io^0, (5) 

where E { is the formation energy of bulk Pu0 2 , namely, E { = |-Ep u o 2 — -^<5-Pu — Eo 2 \- 

To determine reasonable ranges of A/z , the 5-Pu is considered as reference system to 
calculate the formation energy Ef of bulk PuC>2. Since the spin-orbit coupling (SOC) is 
important for certain properties of heavy-metal compounds, we also include SOC effect in 
the calculations of -Ep u o 2 and E S -p u . Finally, we restrict A/x to —4.89 eV ^ A/i ^ 
based on the GGA+U and -4.83 eV ^ A^o < based on the GGA+f/+SOC. The effect 
of spin polarization has been included in calculating Eq 2 . 

III. RESULTS AND DISCUSSION 

A. Surface energy and structural relaxation 

First, the relative stability of the low- index Pu02 surfaces is studied based on the sys- 
tematic calculation of surface energy E s and the detailed analysis of structural relaxation. 
Furthermore, the effects of slab thickness on both surface stability and relaxation are con- 
sidered and discussed. In the following text, we first present the results of non-polar (111) 
and (110) surfaces, and then the polar (001) surface. 

The calculated surface energy for fully relaxed (111) and (110) slabs as a function of the 
thickness is illustrated in Fig. 2(a). Obviously, both the GGA+?7 and LDA+£7 calculations 
give the consistent results that the (110) surface energy is much higher than the (111) 
surface energy. Generally, the calculated E s for (110) is 42% with GGA+U (or 33% with 
LDA+U) higher than for (111). Therefore, the Pu0 2 (lll) surface is much more stable than 
the "more open" (110) surface with relatively higher concentration of the surface dangling 
bonds. For these two non-polar surfaces, Fig. 2(a) furthermore shows two points worthy of 
special notice and further discussion: (i) Despite the large difference in their respective E s 
between GGA+U and LDA+U calculations, the upper value (i.e., the LDA+U result) of the 
(111) surface is notably smaller than the lower value (i.e., the GGA+U result) of the (110) 
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FIG. 2: (Color online) The surface energy E s as a function of the PuC>2 slab thickness: (a) I?g elax 
(solid sphere) of (111) surface, E™ (solid square) and £ , ^ nrelax (open square) of the (110) surface; 
(b) El elax (solid sphere) and E^ nrclaiX (open circle) of defect-(OOl) surfaces. Note that the slab 
thickness is defined by the number N of the Pu-cation layers. Four different terminations of 
defect-(OOl) slabs in (b) are described in Fig. 1(d). 

surface; (ii) The E s of the (111) surface is insensitive to the thickness with steady values 
(0.045eV/A 2 for GGA+U and 0.065eV/A 2 for LDA+U), whereas for the (110) surface the 
evolution of E s as a function of slab thickness shows an oscillating behavior, which indicates 
excellent agreement between LDA+U and GGA+U. 

For the DFT energetic studies of solid materials, it is well known that the GGA calculation 
usually underestimates the experimental value and on the contrary the LDA often reports 
overestimated results for many physical quantities, amongst which the surface energy is 
a typical one. These opposite deviations from the experimental measurement have been 
attributed to the "overbinding" character of LDA and the consequent overcorrection of 



this defect in GGA 



41 



421 ] . However, we are more interested in comparing the relative 



stabilities of different surfaces than in assessing the different performances of LDA and GGA 
functionals, especially in the absence of the experimental data. As far as we are aware, a few 
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existing DFT calculations have given a similar trend in LDA- and GGA-E S results of metal 
oxides. For example, the LDA- E R results are 35% and 22% higher than the GGA results for 
CeO"2(lll) and (110) surfaces 451 ]. respectively, while for MgO(OOl) surface the LDA-.E, is 
25.8% higher than the GGA-E S with the experimental values positioned in between. Here, 
based on the above-mentioned point (i), we are positive that the experimental observation 
will agree on the prominent stability of the PuO^lll) surface. 

We now turn our attention to the evolution with the slab thickness of the surface energies 
in Fig. 2(a). As is known, for a certain cleaved surface, the consequent structural relaxation 
is an effective way to minimize the surface cleavage energy, corresponding to the ^ nrelax ; 
with the contribution defined as the surface relaxation energy AE S =— (El elax — ^ mrelax ). For 
(111) surface, the _E" nrelax (not plotted here) is very close to the E£ elax , producing a quite 
small AE S less than 1.0 meV/A 2 , and the thickness effect on the surface relaxation can be 
neglected. From Fig. 2(a), one can see that the £'^ nrelax of (110) surface is clearly larger 
than the corresponding _£^ clax with their difference (AE S ) being higher than 10 meV/A 2 . 
Additionally, we find that the £^ nrelax is to some extent insensitive to the slab thickness, 
thus the oscillating behavior of ££ elax is tied up to the dependence of AE S upon the slab 
thickness. 

In order to draw a clear comparison of the surface relaxation between (111) and (110) 
surfaces and gain a detailed understanding of the oscillating behavior of (110)-i?g elax , it 
proves to be quite necessary to discuss the surface structures undergoing full relaxations. 
Figures 3(a) and 3(b) show the interlayer relaxations along the directions perpendicular to 
the (111) and (110) surfaces respectively. Here the interlayer relaxation Ai+i^ is given by 
the optimized interlayer distance d- l+ i^ in a relaxed slab compared to the bulk interlayer 
distance d® +1 } along the corresponding direction. Obviously, the signs + and — of Ai+i^ 
indicate expansion and contraction of the interlayer spacing respectively. The stacking 
sequence of the (111) slab with iV=9 consists of 9 "O-Pu-O" blocks, and therefore this (111) 
slab contains totally 27 atomic monolayers (MLs). One can see from Fig. 3(a) that the 
interlayer relaxations are really small, so that the shrinkage ratio of the thickness is only 
~0.25%. The (110) slabs used in Fig. 3(b) contain 9 atomic MLs with two oxygen and 
one plutonium atoms per lxl cell being coplanar. Figure 3(b) shows that (i) the interlayer 
relaxations in the (110) surface region of a few atomic layers are prominently larger than 
those of the (111) surface in Fig. 3(a); (ii) the interlayer relaxations of Pu sublattice are 



10 





0.2r 




0.1 


+ 


0.0 

1 

-0 1 - 



0.1 
0.0 
-0.1 
-0.2 b' 




7\ 



PuOfllO): N=9 



--- A , -O 

i+l,i 

A , -Pu 

i+l,i 

■ i ■ i ■ i . i 



(c) 

surface 



O 



Pll-^Pu-O 



Puyd 

1 Pu-O 



O 



subsurf. 



1 2 3 4 5 6 7 8 
i 

— i — i — i — i — i — i — i — i — i — i — 



0.15 ■ 



9 o.io- 



3 
0. 



0.05 



(d) 



PuO 2 (110): N=3~10 




■ 0.016 <L 

65 

M 

w 

0.014 g 

no 

0.012 1? 



8 9 



10 



FIG. 3: (Color online) The interlayer relaxation Ai+i,; of (a) the (111) slab and (b) (110) slab 
with N = 9. (c) Sketch map of th.6 vertical displacement c?Pu- o in the surface and subsurface of 
the relaxed (110) slab due to the mismatch A2,i-Pu 7^ A2 i i-0 . (d) The dp u -o m both surface and 
subsurface of (110) slab and the relaxation surface energy AE S as a function of the slab thickness. 

larger than the oxygen sublattice, which is especially apparent in the outmost two layers. 
Such mismatch of the relaxations between O- and Pu-sublattices gives rise to the vertical 
displacement <ip u _o between O and Pu atoms which are coplanar in the unrelaxed (110) slab. 
Figure 3(c) gives a sketch map of the <ip u _o in the surface and subsurface layers as a result 
of the mentioned mismatch: A2,i-Pu 7^ A2,i-0. One can see that due to the nonzero dp u _o, 
cation-anion dipoles with inverse orientations are generated in the surface and subsurface 
respectively. Strictly speaking, the relaxed (110) surface is now not a nonpolar surface for 
Pu and O species. As we are aware, this observable surface polarization of PuC^llO) slab 
induced by the structural relaxation is in good agreement with previous DFT calculation 



27| . Besides the interlayer (vertical) relaxation, the inplane (lateral) relaxation of the 
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(110) defect-(OOl) 




FIG. 4: (Color online) The top view of the surface structures of the relaxed (a) (110) slab, and 
the defect-(OOl) slabs: (b) surf-l/surf-2 slabs and (c) surf-l/surf-2 slabs. The white arrows (not 
to scale) indicate the directions of the inpane lateral movements of the surface O atoms in (a) and 
of the subsurface Pu-atoms in (b). 

surface layer (see Fig. 4(a)) tends to shorten the Pu-0 bond on the surface by driving 
two nearest-neighbor oxygen atoms (bonding to the same Pu atom) to close up by ~0.22A, 
leading to the formation of 0-0 dimers on the (110) surface. Furthermore, we have found 
that the vertical displacement dp u _o shows an oscillating behavior with the increasing slab 
thickness while the structure of the inplane 0-0 dimers is to some extent insensitive to the 
slab thickness. To reveal a clear-cut relationship between the structural relaxation and the 
corresponding released energy AE S , we plot the d Pu _ (in both surface and subsurface) and 
AE S as functions of the slab thickness N in Fig. 3(d). One can see that the oscillating 
behaviors of (ip u _o and AE S are quite in-phase, indicating that the oscillations of E s of (110) 
surface as a function of slab thickness originate from the interlayer relaxation. 

The ideal PuO2(001) surface (see Fig. 1(c)) is an unstable polar surface with an over- 
all dipole field. However, it is found that 50% surface O vacancies in our defective (001) 
slab models (see Fig. 1(d)) can effectively quench the dipole field and stabilize the surface. 
Considering that a half oxygen vacancies can usually induce the significant surface recon- 
struction, here we first carry out the first-principles molecular dynamic (FPMD) simulations 
based on GGA+U within the micro-canonical ensemble to sufficiently optimize the defect- 
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(001) surface structures and then calculate their zero-temperature surface energies. From 
the E B result (including _E^ clax and E^ nvelax ) presented in Fig. 2(b) as a function of the slab 
thickness N, one can see that with increasing N, the El elax for the four defect-(OOl) surface 
models converges to ~0.11eV/A 2 , which is still notably larger than that of the (110) surface. 
Interestingly, according to the results of El elax , surf-3 and surf-4 slab models are a bit less 
stable than surf-1 and surf-2 in the whole range of slab thickness that we considered. How- 
ever, the values of ^ nrelax for surf-3 and surf-4 models are lower than those for surf-1 and 
surf-2 models mainly due to the two different distributions of the surface oxygen vacancies, 
namely, the missing-row and uniform types in surf-l/surf-2 and surf-3/surf-4 respectively. 

After the structural optimization by the FPMD simulations, for all four defect- (001) 
surface models with N = 8, the surface oxygens as well as the subsurface oxygens beneath 
relax inward by 0.26 A and 0.14 A respectively, while the subsurface oxygens without surface 
oxygen above relax outward by ~0.17 A. For surf-l/surf-2, the Pu-sublattice relaxes inward 
by ~0.02 A, on the contrary, the Pu-sublattice relaxes outwards by 0.05 A for surf-3/surf- 
4. Accompanied with a slight discrepancy in vertical relaxations of the Pu-sublattice for 
surf-l/surf-2 and surf-3/surf-4, it is found that the difference in AE S for these two types 
of (001) terminations is mainly caused by the distinguishing inplane (lateral) relaxations of 
subsurface Pu-sublattices, which are sketched in Figs. 4(b) and 4(c). For surf-l/surf-2, the 
Pu atoms bonding to the same surface oxygen atom relax to close up by ~0.45 A and this 
periodic lattice distortion consequently provoke the zigzag manner reconstruction of surface 
oxygen-lattice from the initial linear chain. Interestingly, this peculiar reconstruction was 
experimentaly observed in the defective polar (001) surface of uranium dioxide UO2 with 
50% oxygen vacancies [4^. However, because of the uniform distribution of the surface 
O- vacancies, the Pu and O atoms in surf-3/surf-4 keep lateral inaction. For the defective 
(001) surfaces, our current results show that both the slab thickness and the distributions of 
the surface O-vacancies can notably impact the surface stability, and there may be several 
more stabilizing mechanisms coexisting on the polar (001) surface. 

To briefly summarize our results in this section, we give the relative order of stability for 
low-index PuC>2 surfaces, namely, (111) > (110 ) > defect-(OOl) > polar-(OOl), which is well 



consistent with that of CeC>2 
as PuC>2. 
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45) and U0 2 



46j , which are of the same fluorite structure 
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B. Surface electronic structure and work function 



The bulk P11O2 is considered to be the AFM Mott insulator according to the experimental 



report 



121 ] . Here the atom-projected density of the electronic states (PDOSs) for the Pu and 



oxygen atoms in the bulk PuC>2 and on the relaxed (111), (HO), and defect-(OOl) surfaces are 
shown in Fig. 5. Since the GGA+U and LDA+£7 give the similar description of the PDOS, 
here we only plot the GGA+U results. The orbital-resolved PDOS of the bulk PuC>2 at the 



ground state has been ca.cu.ated and ana. y zed in detai. by p r ev,ous DFT+f/ UMMM 
and hybrid DFT [9|, |16|] studies, and those theoretical results of DOS are usually tested by 
comparing with the experimental photoelectron spectroscopy (PES) measurements 

aa. 

Here, we replot the PDOS of the bulk PUO2 with AFM phase as a benchmark for those 
of the Pu-0 atoms on different Pu02 surfaces, aiming at finding significant influence in the 
electronic structure by the inclusion of surface effect. The PDOS of bulk PUO2 in Fig. 5(a) 
demonstrates the following features: (i) Above the Fermi level, the insulatingband gap is 
about 1.8 eV, which is in good agreement with the experimental measurement [6(; (ii) Below 
the Fermi level the highest occupied band (HOB) with a range of —5 to eV is mainly 
the 5/(Pu)-2p(0) hybridization, with little contributions from 6p and 6d states of Pu; (iii) 
The lower occupied band (LOB) with a range of —21 to —13 eV consists of Pu-6p and 0-2s 
states. 

From the surface DOS in Fig. 5, one can see that the PDOS distribution for the (111) 
surface shows a close resemblance to the case of the bulk Pu02, specifically, the similarities 
in the insulating band gap and the structures of Pu-5/ state with two pronounced peaks are 
so strong that the slight contraction and shift-up of the HOB are almost covered up. For the 
(110) surface, the insulating band gap reduces to ~1.6 eV, and particularly, the two-peak 
structure of Pu-5/ disappears mainly due to the existing surface polarization with nonzero 
^Pu-o? which modifies the crystal symmetry of the oxide surface layer. For the (001) surface, 
since the surface layer of surf-1 slab model used in Fig. 5(d) is not the stoichiometric Pu02 
but the 'Pu-O', the insulating band gap further reduces to ~1.4eV, and a sharp peak of Pu- 
5/ state emerges below the Fermi level, which implies the increase in the localized correlation 
of the Pu-5/ electronic state due to the presence of oxygen vacancies. These facts suggest 
that the surface effect of Pu02 upon the electronic structure of the bulk phase appears to 
be insignificant for the stable (111) surface, to some extent significant for the (110) surface, 
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FIG. 5: (Color online) The atom-projected (orbital resolved) DOS for (a) bulk Pu02, (b) (111) 
surface, (c) (110) surface, and (d) defect-(OOl) surface. The Fermi level is indicated by the vertical 
dashed line at eV. 
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FIG. 6: (Color online) The calculated work function $ as a function of the slab thickness: (a) 
PuO"2(lll) and (110) surfaces, (b) the ideal (001) and the four defect-(OOl) surfaces. 

and remarkable for the defect-(OOl) surface. In order to be able to theoretically reproduce 
the whole PES spectra of PuC>2 layers by the right description of the complex behavior of 
Pu-5/ state, there is clearly much left to be done. 

In addition to the PDOS, we have also calculated the work function $ of low-index PuC>2 
surfaces, and plotted them in Fig. 6. The work function $ is the minimum energy required to 
remove an electron from the surface to the vacuum and can be written as $ = Vacuum - E-p, 
where Kacuum is the planar-averaged electrostatic potential in the middle of the vacuum and 
Ep is the Fermi energy of the system. Therefore, work function is one fundamental physical 
quantity for the surface reactivity. Furthermore, a modified or tunable $ can be useful for 
applications such as catalysis, because a slight change in the energy scale is exponentially 
amplified for chemical reactions. 

From Fig. 6(a), one can see that the <3> of stable (111) surface with an average value of 
~6.1 eV (GGA+L7) or ~6.3 eV (LDA+U) is much higher than that of the (110) surface with 
an average value of ~4.7 eV (GGA+C/) or ~4.8 eV (LDA+U) . Thus, for nonpolar surfaces, 
stable (111) surface will show its inertness in the surface chemical reactions, and the more 
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open (110) surface is expected to be chemically active. Figure 6(a) also demonstrates the 
convergence behavior of $ as a function of the slab thickness. Here, it is found that the 
GGA+U and LDA+[7 results of $ are in general agreement. For the stable (111) surface, 
the work function shows less responsiveness to the thickness effect, as well as its surface 
energy E s in Fig. 2(a). On the contrary, both $ and E s of the (110) surface are sensitive 
to the slab thickness with a convergent oscillation. Combining with the surface relaxation 
results given in Fig. 3(d), one can conclude that the thickness effect modifies both the 
surface stability and surface chemical activity through the structural relaxations. Usually, 
the thickness of the oxide film formed on Pu metal during stockpile process is typically of 
nanometer scale. Thus, our present finding of thickness-selective surface activity may help to 
deepen the understanding of the microscopic mechanisms for the chemical reaction of small 
molecules (such as water) on oxidized Pu surfaces, which is fundamental to the safety issue 
of nuclear industry. Interestingly, a recent DFT study has reported that the thickness 
effect of MgO film can be used to control the dissociation of water molecule on surface. 

For the case of the polar (001) surface shown in Fig. 5(b), the work function is mainly 
dominated by the strength of the anion-cation dipole, which impedes the escape of electrons. 
Therefore, one can see that the work function of the ideal (001) surface is close to 8 eV, which 
is much higher than that of the defect-(OOl) surfaces with an average value of ~6.4 eV. Due 
to the significant influence of the existing dipole on the (001) surface, it is not reasonable to 
compare its surface chemical activity with ideal (110) or (111) non-polar surface merely via 
the surface work function. 



C. Effect of oxygen vacancy 

In this section, we focus on the effect of O- vacancy with various concentrations upon the 
surface activity and surface relative stability by using the static GGA+U calculation and 
the approach of u ab initio atomistic thermodynamics". Here, our current study is mainly 
driven by the following motivations: (i) to explain the difference in the surface chemical 
activity between PuC>2 and a-Pu203, as mentioned in Sec-I; (ii) to discuss the mechanism 
of creating surface oxygen- vacancy in the cancelation of the polarity for PuC^OOl) surface; 
and (iii) to explore the stable surface phase of Pu02 in an oxidizing environment. Amongst 
these listed issues, the O-vacancy is obviously the major factor. 
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TABLE I: The calculated work function $ (in eV) of Pu0 2 (lll), (110) and (001) surfaces. 







C v =0 


C v =l/9 


C v =l/4 


C v =l/2 


C v =3/4 


C v =l 


(Ill): 


on-surface 


6.19 


5.09 


4.35 


4.07 






(111): 


subsurface 


6.19 


5.36 


5.18 


4.87 


4.49 


2.57 


(110): 


on-surface 


4.70 




3.84 


3.57 


2.80 


2.44 


(001): 


on-surface 


7.93 




7.24 


6.63 


4.68 


3.25 



In the calculation, to eliminate the thickness effect of the Pu0 2 slab, we employ (111), 

(110) , and (001) slabs with N — 6, 8, and 8 respectively. The various concentrations of 
oxygen vacancy (Cy) are modeled by removing different amounts of oxygen atoms from 
ideal slabs with a series of surface unit cells. Here, Cy is the ratio between the number of 
O-vacancies and the total number of O atoms on the ideal surface layer. Specifically, for the 

(111) surface, slabs of (2 x 2) and (3 x 3) unit cells are employed to create Cy of |, |, |, |, 
|, and 1.0. For the (110) surface, slab of (1 x 2) unit cells is used to create Cy of |, |, |, 
and 1.0. For the (001) surface, slab of (y/2 x y/2)RA5° unit cells can create Cy of \, |, §, 
and 1.0. 

Before the discussion on the evolution of work function as a function of Cy, we first present 
the O-vacancy effect upon the structural relaxation of the (111) slab. The static calculation 
demonstrates that the (111) surface with on-surface O-vacancy is slightly preferred in total 
energy when Cy ^1/2. However, when 1/2<CV ^1, the subsurface oxygen atoms, each 
of which sharing the same unit surface cell with one certain surface oxygen, will break 
through the above Pu-terminated layer to form a complete O-terminated surface, leaving 
the subsurface O-vacancies at their former sites. Therefore, when C v >l/2, all the on-surface 
O-vacancies will convert to be the subsurface O-vacancies by a significant reconstruction. 

The GGA+£7 calculated work function $ of (111), (HO), and (001) surfaces for different 
values of Cy is listed in Table I, where the "on-surface" and the "subsurface" denote the 
initial pure on-surface and subsurface distributions of O-vacancies for the (111) surface. One 
can see from Table I that for all three surfaces the work function will monotonically reduce 
with increasing Cv. Therefore, the introduction of O-vacancy can prominently enhance 
the surface chemical activity of non-polar (111) and (110) surfaces. For instance, the work 
function of the ideal (111) surface is 6.19 eV, while a low Cy = 1/9 (1/4) of on-surface 
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FIG. 7: (Color online) Surface free energy difference A7 of (a) Pu0 2 (lll), (b) (110) and (c) (001) 
surfaces with various concentrations of O-vacancy Cv as a function of the oxygen chemical potential 
A/xo; with the corresponding pressure lines at T=300 K and T=600 K. The low-energy surface 
terminations are drawn with the thick lines in (a)-(c), which are gathered in (d) with the ideal 
(111) surface as the reference configuration. 

O-vacancy can effectively depress the work function by 1.1 (1.84) eV and efficiently amplify 
the probability of the chemical reaction between the Pu02(lll) surface and other small 
gaseous molecules such as H2 and H2O, which will be investigated in our next work. This 
result can be also extended to explain the significant difference in the chemical activities 
between PuC>2 and a-Pu203, the latter has been found in experiment to be more active 
in interacting with small molecules. 

Assuming that the PuC>2 surface is in equilibrium with an external oxygen environment 
and translating the oxygen chemical potential into temperature and pressure conditions 
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according to Eq. (3) and Eq. (4) respectively, we first discuss effect of the O-vacancy 
upon the relative stability of one certain surface. Figures 7(a)- (c) present the Gibbs surface 
free energy difference A7 of (111), (HO) and (001) surfaces respectively. Here the A7 is 
calculated by Eq. (3) with the corresponding ideal surface as the reference system. For 
the (HI) surface, from Fig. 7(a) one can see that (i) the ideal, vacancy-free (111) surface 
is the most stable structure under the oxygen rich conditions with A/io ^ —2.49 eV; (ii) 
then defect- (111) surface with on-surface O-vacancy of Cy — 1/9 becomes stable within 
a very narrow range of —2.63 ^ A/io < —2.49 eV; (iii) for further reducing environment 
with Afi ^ —2.63 eV, the defect-(lll) surface with subsurface O-vacancy of high C v =8/9 
becomes the most stable. 

For the (110) surface, Fig. 7(b) shows that (i) the ideal (110) surface is most stable 
within the range of the A/io ^ —1.96 eV; (ii) the defect surface with Cy—l/A becomes most 
stable within —2.54 < A/i < —1.96 eV; (iii) the defect surface with Cy—1/2 in succession 
becomes most stable when —2.93 ^ A/x < —2.54 eV; (iv) the defect surface with Cy— 3/4 
becomes most stable when —3.50 ^ A/i < —2.93 eV, and the Pu-terminated (110) surface 
with Cy— 1 is not a stable surface phase in the whole range of the allowed A/io- 

From the results of polar-(OOl) surface in Fig. 7(c), we have found that (i) under the O- 
rich conditions, the ideal (001) surface is unstable, however the 50% surface O- vacancies can 
efficiently stabilize the polar surface (ii) when A/i Q ^ —1.69 eV, the defect-(OOl) surface with 
CV=3/4 is the optimal surface structure; (iii) within the whole allowed range of A/io, the 
Pu-terminated (001) surface is unstable so as the vacancy-free O-terminated (001) surface. 

We collect all the stable surface phases from Fig. 7(a)- (c), and summarize them in Fig. 
7(d) by taking the ideal (111) surface as the reference structure. One can see that the ideal 
(111) surface, defect-(lll) surface with low on-surface O-vacancy concentration of Cy— 1/9 
and with high sub-surface Cy =8/9 are the stable surface structures. That is to say, the 
ideal (111) surface is stable under the oxygen-rich conditions, while for an oxygen- reducing 
environment the (111) surface with nearly one monolayer subsurface oxygen removed become 
stable, and the on-surface oxygen vacancy with low Cy of 1/9 can minimize the Gibbs surface 
energy in a very narrow range of A/iQ. 

In Table II, we have listed the O-vacancy formation energies Ey, which can be defined as 
Ev = Wv~ [-^s d iab ct ~ E lub + ^o-v • \E 02 ~\ , where iVo-v is the total number of the O-vacancy 
in a defective slab, E^ ct , E 1 ^ 1 and E Q2 are the total energies of the defective slab, ideal 
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TABLE II: The calculated formation energies of the O-vacancy E v (in eV). 
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1.96 


2.54 


2.93 




3.50 
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on-surface 




-3.38 


-1.93 


-0.73 




1.11 



slab and a free oxygen molecule, respectively. One can see for the (111) surface that the Ey 
does not show considerable change except in the extreme case of Cy = 1 with a maximum 
value Ey— 3.20 eV. On the contrary, the Ey of the (110) surface is sensitive to the Cy, 
namely, the Ey monotonically increases with increasing Cy, indicating a notable interaction 
between the vacancies. Finally, the polar (001) surface is a special case. The minus Ey 
indicates that the formation of surface O-vacancy is an exothermic process, and at the same 
time stabilizes the polar surface. However the Ey also monotonically increases with the Cy 
and rises to +1.11 eV when Cy = 1. 

IV. CONCLUSIONS 

To conclude, we have systematically studied the basic surface properties of low-index 
Pu02(lll), (HO), and (001) surface by means of the first-principles DFT calculations within 
the LDA+£7 and GGA+U frameworks. The defect-free O-terminated (111) surface is found 
to be most stable, possessing the lowest E s that is insensitive to the thickness of the film. 
The surface energy of the non-polar (110) surface is 33% to 42% higher than that of the 
(111) surface, accompanying with an oscillating behavior with the film thickness. The polar 
(001) surface has been modeled using 50% oxygen vacancies to cancel the polarity. The 
residual surface oxygen atoms have been found to reconstruct in a zigzag manner along the 
<100> direction. In connected with the relative order of stability for these three low-index 
surfaces, our calculated surface electronic structures have displayed from insignificant to 
remarkable deviation from the bulk case. The work function $ has also been systematically 
investigated, and a high value of about 6.19 eV for the most stable (111) surface indicates its 
low chemical activity. Remarkably, this value can be reduced to 4.35 eV with 25% oxygen- 
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vacancy present on the surface. This conclusion can be used to explain the difference in the 
surface chemical activities between PuC>2 and a-Pu2C>3. 

We have also investigated the surface thermodynamics in an oxygen environment. Our 
results have indicated that under oxygen-rich conditions, the stoichiometric (111) surface 
is most stable. Under oxygen-reducing conditions, the on-surface O-vacancy of low concen- 
tration Cy = 1/9 can slightly minimize the Gibbs surface energy 7 of (111) in a narrow 
range of the oxygen chemical potential A/xo- For the highly reducing conditions, the (111) 
surface with nearly one monolayer subsurface oxygen removed (Cy = 8/9) becomes most 
stable, where the upper layers resemble the /3-Pu 2 O 3 (0001) surface. Based on these system- 
atic results, our current study may provide a guiding line to understand various chemical 
properties and processes occurred on Pu0 2 surfaces. 
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